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Abstract 

The AGBNP2 implicit solvent model, an evolution of the Analytical Generalized 
Bom plus Non-Polar (AGBNP) model we have previously reported, is presented with 
the aim of modeling hydration effects beyond those described by conventional con- 
tinuum dielectric representations. A new empirical hydration free energy component 
based on a procedure to locate and score hydration sites on the solute surface is in- 
troduced to model first solvation shell effects, such as hydrogen bonding, which are 
poorly described by continuum dielectric models. This new component is added to 
the Generalized Born and non-polar AGBNP models which have been improved with 
respect to the description of the solute volume description. We have introduced an 
analytical Solvent Excluding Volume (SEV) model which reduces the effect of spurious 
high-dielectric interstitial spaces present in conventional van der Waals representations 
of the solute volume. The AGBNP2 model is parametrized and tested with respect 
to experimental hydration free energies of small molecules and the results of explicit 
solvent simulations. Modeling the granularity of water is one of the main design princi- 
ples employed for the the first shell solvation function and the SEV model, by requiring 
that water locations have a minimum available volume based on the size of a water 
molecule. We show that the new volumetric model produces Born radii and surface 
areas in good agreement with accurate numerical evaluations. The results of Molecular 
Dynamics simulations of a series of mini-proteins show that the new model produces 
conformational ensembles in substantially better agreement with reference explicit sol- 
vent ensembles than the original AGBNP model with respect to both structural and 
energetics measures. 



1 Introduction 

Water plays a fundamental role in virtually all biological processes. The accurate modeling 
of hydration thermodynamics is therefore essential for studying protein conformational equi- 
libria, aggregation, and binding. Explicit solvent models arguably provide the most detailed 
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and complete description of hydration phenomena.- They are, however, computationally de- 
manding not only because of the large number of solvent atoms involved, but also because of 
the need to average over many solvent configurations to obtain meaningful thermodynamic 
data. Implicit solvent models,- which are based on the statistical mechanics concept of the 
solvent potential of mean force,- have been shown to be useful alternatives to explicit solva- 
tion for applications including protein folding and binding,- and small molecule hydration 
free energy prediction.- 

Modern implicit solvent modelsSiZ include distinct estimators for the non-polar and elec- 
trostatic components of the hydration free energy. The non-polar component corresponds 
to the free energy of hydration of the uncharged solute while the electrostatic component 
corresponds to the free energy of turning on the solute partial charges. The latter is typically 
modeled treating the water solvent as a uniform high-dielectric continuum.- Methods based 
on the numerical solution of the Poisson-Boltzmann (PB) equation- provide a virtually exact 
representation of the response of the solvent within the dielectric continuum approximation. 
Recent advances extending dielectric continuum approaches have focused on the develop- 
ment of Generalized Born (GB) models,-- which have been shown to reproduce with good 
accuracy PB and explicit solvent-i^i results at a fraction of the computational expense. The 
development of computationally efficient analytical and different iable GB methods based 
on pairwise descreening schemes-'^'^ has made possible the integration of GB models in 
molecular dynamics packages for biological simulations.—'^'^ 

The non-polar hydration free energy component accounts for all non-electrostatic solute- 
solvent interactions as well as hydrophobic interactions,— which are essential driving forces 
in biological processes such as protein folding— li^i^i^i and binding.— 12^'2^'2^ Historically the 
non-polar hydration free energy has been modeled by empirical surface area models^ which 
are still widely employed.— i^IiS^iSiiSSiSiiS^i^i^i^ Surface area models are useful as a first ap- 
proximation, however qualitative deficiencies have been observed.— '^'^'^'^'^i^ 

Few years ago we presented the Analytical Generalized Born plus Non-Polar (AGBNP) 
implicit solvent model, --^ which introduced two key innovations with respect to both the 
electrostatic and non-polar components. Unlike most implicit solvent models, the AGBNP 
non-polar hydration free energy model includes distinct estimators for the solute-solvent van 
der Waals dispersion energy and cavity formation work components. The main advantages of 
a model based on the cavity/ dispersion decomposition of the non-polar solvation free energy 
stem from its ability to describe both medium range solute-solvent dispersion interactions, 
which depend on solute composition, as well as conformational equilibria dominated by 
short-range hydrophobic interactions, which mainly depend on solute conformation alone.— 
A series of studies highlight the importance of the balance between hydrophobicity and 
dispersion interactions in regulating the structure of the hydration shell and the strength 
of interactions between macromolecules.— i^i^ In AGBNP the work of cavity formation is 
described by a surface area-dependent model,— '^'^'^ while the dispersion estimator is based 
on the integral of van der Waals solute-solvent interactions over the solvent modeled as a 
uniform continuum.— This form of the non-polar estimator had been motivated by a series of 
earlier studies^i^i^i^i^i^ and has since been shown by u^sSiSs^iSS others^'^'^'^ to be 
qualitatively superior to models based only on the surface area in reproducing explicit solvent 
results as well as rationalizing structural and thermodynamical experimental observations. 

The electrostatic solvation model in AGBNP is based on the pairwise descreening GB 
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scheme^ whereby the Born radius of each atom is obtained by summing an appropriate de- 
screening function over its neighbors. The main distinction between the AGBNP GB model 
and conventional pairwise descreening implementations is that in AGBNP the volume scaling 
factors, which offset the overcounting of regions of space occupied by more than one atom, 
are computed from the geometry of the molecule rather than being introduced as geometry- 
independent parameters fit to either experiments or to numerical Poisson-Boltzmann re- 
sults.— i^i^i^ The reduction of number of parameters achieved with this strategy improves 
the transferability of the model to unusual functional groups often found in ligand molecules, 
which would otherwise require the introduction of numerous parameters.— 

Given its characteristics, the AGBNP model has been mainly targeted to applications in- 
volving molecular dynamics canonical conformational sampling, and to the study of protein- 
ligand complexes. Since its inception the model has been employed in the investigation 
of a wide variety of biomolecular problems ranging from peptide conformational propen- 
sity prediction and folding,^'^'^!^ ensemble-based interpretation of NMR experiments,—!^ 
protein loop homology modeling,— the study of intrinsically disordered proteins,^^ ligand- 
induced conformational changes in proteins,—*^ conformational equilibria of protein-ligand 
complexes,—!^ protein-ligand binding affinity prediction,— and structure-based vaccine de- 
sign.— The AGBNP model has been re-implemented and adopted with minor modifications 
by other investigators.-^'^ The main elements of the AGBNP non-polar and electrostatic 
models have been independently validated,— "^"^"^ and have been incorporated in recently 
proposed hydration free energy models.—"^ 

In this work we present a new implicit solvent model named AGBNP2 which builds upon 
the original AGBNP implementation (hereafter referred to as AGBNP 1) and improves it 
with respect to the the description of the solute volume and the treatment of short-range 
solute-water electrostatic interactions. 

Continuum dielectric models assume that the solvent can be described by a linear and 
uniform dielectric medium.— This assumption is generally valid at the macroscopic level, 
however at the molecular level water exhibits significant deviations from this behavior.- 
Non-linear dielectric responses, the non-uniform distribution of water molecules, charge 
asymmetry and electrostriction effect&Sfi are all phenomena originating from the finite size 
and internal structure of water molecules as well as their specific interactions which are 
not taken into account by continuum dielectric models. Some of these effects are qualita- 
tively captured by standard classical fixed-charge explicit water models, however others, such 
as polarization and hydrogen bonding interactions, can be fully modeled only by adopting 
more complex physical models.— GB models make further simplifications in addition to the 
dielectric continuum approximation, thereby compounding the challenge of achieving with 
GB-based implicit solvent models the level of realism required to reliably model phenomena, 
such as protein folding and binding, characterized by relatively small free energy changes. 

In the face of these challenges a reasonable approach is to adopt an empirical hydra- 
tion free energy model motivated by physical arguments^^ parametrized with respect to 
experimental data on small molecule hydration.— Models of this kind typically score con- 
formations on the basis of the degree of solvent exposure of solute atoms. Historically— the 
solvent accessible surface area of the solute has been used for this purpose, while modern 
implementations suitable for conformational sampling applications often employ computa- 
tionally convenient volumetric measures.— In this work we take this general approach but 
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we retain the GB model component which we assume to form a useful baseline to describe 
the long-range influence of the water medium. The empirical parametrized component of the 
model then takes the form of an empirical first solvation shell correction function designed 
so as to absorb hydration effects not accurately described by the GB model. Specifically, as 
described below, we employ a short-range analytical hydrogen bonding correction function 
based on the degree of water occupancy taking into account the finite size of water molecules, 
of appropriately chosen hydration sites on the solute surface. The aim of this model is to ef- 
fectively introduce some explicit solvation features without actually adding water molecules 
to the system as for example done in hybrid explicit /implicit models.—*^ 

In this work we also enhance the volume description of the solute, which in AGBNPl is 
modeled by means of atomic spheres of radius equal to the atomic van der Waals radius. The 
deficiencies of the van der Waals solute volume model have been recognized.— They stem 
from the presence of high-dielectric interstitial spaces in the solute interior which are too 
small to contain discrete water molecules. These spurious high dielectric spaces contribute 
to the hydration of buried or partially buried atoms causing underestimation of desolvation 
effects. The volume enclosed by the molecular surface (MS), defined as the surface produced 
by a solvent spherical probe (of 1.4 A radius for water) rolling on the van der Waals surface 
of the solute,-^ represents the region which is inaccessible to water molecules and is often 
referred as the Solvent Excluding Volume (SEV).-2S The SEV, lacking the spurious high- 
dielectric interstitial spaces, provides a better representation of the low-dielectric region 
associated with the solute. For this reason accurate Poisson-Boltzmann solver s,-'^'^^ have 
employed the SEV description of the solute region. 

Despite its clear advantages, the lack of of analytical and computational efficient de- 
scriptions of the SEV have hampered its deployment in conjunction with GB models for 
molecular dynamics applications. The GBMV series of models^>2^i2i achieve high accuracy 
relative to numerical Poisson calculations in part by employing the SEV description of the 
solute volume. The analytical versions of GBMV—^ describe the SEV volume by means 
of a continuous and differentiable solute density function which is integrated on a grid to 
yield atomic Born radii. In this work we present a model for the SEV that preserves the 
analytical pairwise atomic descreening approach employed in the AGBNPl model,— which 
avoids computations on a grid. We show that this approximate model reproduces all of the 
key features of the SEV while yielding the same favorable algorithmic scaling of pairwise 
descreening approaches. 

2 Methods 

2.1 The Analytical Generalized Born plus Non-Polar Implicit Sol- 
vent Model (AGBNP) 

In this section we briefly review the formulation of the AGBNPl implicit solvent model, 
which forms the basis for the new AGBNP2 model. A full account can be found in the 
original reference.^ The AGBNPl hydration free energy AG^ is defined as 
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— AGelec + AGcav + AGvdW, (1) 

where AGeiec is the electrostatic contribution to the solvation free energy and AGnp includes 
non-electrostatic contributions. AGnp is further decomposed into a cavity hydration free 
energy AGcav and a solute-solvent van der Waals dispersion interaction component AGvdw- 

2.1.1 Geometrical estimators 

Each free energy component in Eq. ([1]) is ultimately based on a analytical geometrical de- 
scription of the solute according to which the solute volume is modeled as a set of overlapping 
atomic spheres of radii Ri centered on the atomic positions r,. Hydrogen atoms do not con- 
tribute to the solute volume. The solute volume is computed using the Poincare formula (also 
known as the inclusion-exclusion formula) for the volume of the union of a set of intersecting 
elements 

^ = E^^-E^^.+ E v,^,-... (2) 

i i<j i<j<k 

where Vi = AttR'^/S is the volume of atom i, Vij is the volume of intersection of atoms i and 
j (second order intersection), Vijk is the volume of intersection of atoms i, j, and k (third 
order intersection), and so on. The overlap volumes are approximated by the overlap integral, 
VjI n5 available in analytical form, between n Gaussian density functions each corresponding 
to a solute atom: 



12. ..n 



J dhp,{v)p,{v)...p^{v), (3) 
where the Gaussian density function for atom i is 

Pi(r)=pexp -Ci{r-rif , (4) 



where 



R, 



(5) 
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and 

and K is a dimensionless parameter that regulates the diffuseness of the atomic Gaussian 
function. In the AGBNPl formulation k = 2.227. 

Since Gaussian integrals are in principle non-zero for any finite distances between the 
Gaussian densities, the question arises of how to implement in practice Eq. ([2]) without 
having to compute the combinatorial large number of mostly small overlap volumes between 
all of the atoms of the molecule. Although not mentioned in reference,^ this problem had 
been addressed in AGBNPl by introducing a switching function that smoothly reduces to 
zero the overlap volume between two or more Gaussians when the overlap volume is smaller 
than a certain value. If V12... is the value of the Gaussian overlap volume between a set of 
atoms, the corresponding overlap volume V12... used in Eq. ([2]) is set as 

] nl... < vi 

Vi2...n = { Vl2...nU{u) Vi < Vf^.^n < V2 , (7) 
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where 

V&... - vi ..^ 
u = , 8 

/^(x) =x^(10-15x + 6x^), (9) 

and Vi = 0.2 A'^ and f 2 = 2 A^. This effectively sets to zero Gaussian overlap volumes smaller 

than f 1 , leaves volumes above V2 unchanged and smoothly reduces volumes in between these 

two limits. This scheme drastically reduces the number of overlap volumes that need to 

be calculated since the fact that an n-body overlap volume Vi2...n between n atoms is zero 

guarantees that all of the n + 1-body overlap volumes corresponding to the same set of 

atoms plus one additional atom are also zero. (Note below that the formulation of AGBNP2 

employs smaller values of Vi and V2 to improve the accuracy of surface areas). 

The van der Waals surface area Ai of atom i, which is another geometrical descriptor of 

the model, is based on the derivative dV/ dRi of the solute volume with respect to the radius 
p 95 



where V is given by Eq. ([2]) and 



X > 



with a = 5 A^, is a filter function which prevents negative values for the surface areas for 
buried atoms while inducing negligible changes to the surface areas of solvent-exposed atoms. 
The model further defines the so-called self- volume VI of atom i as 

v: = v.-\Y.y^3 + \Y.yr^'^+--- ■ (12) 

which is computed similarly to the solute volume and measures the solute volume that is 
considered to belong exclusively to this atom. Due to the overlaps with other atoms, the 
self-volume VI of an atom is smaller than the van der Waals volume Vi of the atom. The 
ratio 

V 

= ^ < 1 (13) 
is a volume scaling factor is used below in the evaluation of the Born radii. 

2.1.2 Electrostatic Model 

The electrostatic hydration free energy is modeled using a continuous dielectric representa- 
tion of the water solvent using the Generalized Born (GB) approximation 

i i<j Jij 

where 
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and where ein is the dielectric constant of the interior of the solute, is the dielectric 
constant of the solvent; and qj are the charges of atom i and j, and 



= + B,Bj exp(-4/4i?,5,) . (16) 

In Eqs. f ll4l) - (fT6l) Bi denotes the Born radius of atom i which, under the Coulomb field 
approximation, is given by the inverse of the integral over the solvent region of the negative 
4th power of the distance function centered on atom 

P^ = l/B,= ^^ dh^^-^^. (17) 

47r JSolvent (r — Ti)^ 

In the AGBNPl model this integral is approximated by a so-called pairwise descreening 
formula ^ ^ 

where Ri is the van der Waals radius of atom i, sji is the volume scaling factor for atom 
j [Eq- (1131) ] when atom i is removed from the solute and Qji is the integral (available in 
analytic form see Appendix X of reference X) of the function (r — r^)"^ over the volume of 
the sphere corresponding to solute atom j that lies outside the sphere corresponding to atom 
i. Eq. f|T8l) applies to all of the atoms i of the solute (hydrogen atoms and heavy atoms) 
whereas the sum over j includes only heavy atoms. The AGBNPl estimate for the Born 
radii Bi are finally computed from the inverse Born radii f3i from Eq. (ITSl) after processing 
them through the function: 

= m) = { v"^^ > V (19) 

where b^' = 50 A. The filter function Eq ( |T9i) is designed to prevent the occurrence of 
negative Born radii or Born radii larger than 50 A. The goal of the filter function is simply 
to increase the robustness of the algorithm in limiting cases. The filter function has negligible 
effect for the most commonly observed Born radii smaller than 20 A. 

In the AGBNPl model the scaling factors sji are approximated by the expression 

+ (20) 

where Sj is given by Eq. f[T^ and Vij is the two body overlap volume between atoms i and 
j. Also, in the original AGBNP formulation the computation of the scaling factors and the 
descreening function in Eq. f|T8|) employed the van der Waals radii for the atoms of the solute 
and the associated Gaussian densities. These two aspects have been modified in the new 
formulation (AGBNP2) as described below. 

2.1.3 Non-Polar Model 

The non-polar hydration free energy is decomposed into the cavity hydration free energy 
AGcav and the solute-solvent van der Waals dispersion interaction component AGvdw- The 
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cavity component is described by a surface area model— i^i^i^ 

AGcav = E7i^^' (21) 

i 

where the summation runs over the solute heavy atoms, Ai is the van der Waals surface area 
of atom i from Eq. fllUp . and 7^ is the surface tension parameter assigned to atom i (see 
Table X of reference X). Surface areas are computed using augmented radii for the atoms 
of the solute and the associated Gaussian densities. Augmented radii are defined as the the 
van der Waals radii (Table X of reference X) plus a 0.5 A offset. The computation of the 
atomic surface areas in AGBNP2 is mostly unchanged from the original implementation, ^2, 
with the exception of the values of the switching function cutoff parameters vi and V2 of 
Eq. (171), which in the new model are set as vi = 0.01 A'^ and V2 = 0.1 A^. This change was 
deemed necessary to improve the accuracy of the surface areas which in the new model also 
affect the Born radii estimates through Eq. (128!) below. 

The solute-solvent van der Waals free energy term is modeled by the expression 

AG.dw = E«^(5^+\)3' (22) 

where is an adjustable dimensionless parameter on the order of 1 (see Table X of reference 
X) and 

16 

o-i = — —Tipweiwcrf^ , (23) 

where pw = 0.033428 A~'^ is the number density of water at standard conditions, and 
and eiuj are the OPLS force field— Lennard- Jones interaction parameters for the interaction 
of solute atom i with the oxygen atom of the TIP4P water model.— If cxj and are the 
OPLS Lennard- Jones parameters for atom i, 

(7iw = ^/(Ji(7w (24) 

^iw = V ^i^w 5 (25) 

where = 3.15365 A, and = 0.155 kcal/mol are the Lennard- Jones parameters of the 
TIP4P water oxygen. In Eq. fl22l) Bi is the Born radius of atom i from Eqs. flTSl) and flTOl) 
and Rw = 1.4 A is a parameter corresponding to the radius of a water molecule. 



2.2 Pairwise Descreening Model using the Solvent Excluding Vol- 
ume 

When using van der Waals radii to describe the solute volume, small crevices between atoms 
(Fig. m panel A) are incorrectly considered as high-dielectric solvent regions,— "^"^ leading to 
underestimation of the Born radii, particularly for buried atoms. The van der Waals volume 
description implicitly ignores the fact that the finite size of water molecules prevents them 
to occupy sites that, even though are not within solute atoms, are too small to be occupied 
by water molecules. Ideally a model for the Born radii would include in the descreening 
calculation all of the volume excluded from water either because it is occupied by a solute 
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atom or because it is located in an interstitial region inaccessible to water molecules. We 
denote this volume as the Solvent Excluding Volume (SEV). A realistic description of the 
SEV is the volume enclosed within the molecular surface^ of the solute obtained by tracing 
the surface of contact of a sphere with a radius characteristic of a water molecule (typically 
1.4 A) rolling over the van der Waals surface of the solute. The main characteristic of this 
definition of the SEV (see Fig [2]) is that, unlike the van der Waals volume, it lacks small 
interstitial spaces while it closely resembles the van der Waals volume near the solute-solvent 
interface. The molecular surface description of the SEV can not be easily implemented into 
an analytical formulation. In this section we will present an analytical description of the SEV 
for the purpose of the pairwise descreening computation of the Born radii as implemented 
in AGBNP2 that preserves the main characteristics of the molecular surface description of 
the SEV. 

The main ideas underpinning the SEV model presented here are illustrated in Fig. [1] 
We start with the van der Waals representation of the solute (model A) which presents an 
undesirable high dielectric interstitial space between the two groups of atoms. Increasing the 
atomic radii leads to a representation (model B) in which the interstitial space is removed 
but that also incorrectly excludes solvent from the surface of solvent exposed atoms. This 
representation is therefore replaced with one in which the effective volume of each atom in 
B is reduced by the volume subtended between the solvent-exposed surface of each atom 
and its van der Waals radius (panel C). This process yields model D in which the effective 
volume of the most buried atom is larger than those of the solvent-exposed atoms. This SEV 
model covers the interstitial high-dielectric spaces present in a van der Waals description 
of the solute volume, while approximately maintaining the correct van der Waals volume 
description of atoms at the solute surface as in the molecular surface description of the SEV 
(Fig. [2D. 

These ideas have been implemented in the AGBNP2 model as follows. The main modi- 
fication consists of adopting for the pairwise descreening Generalized Born formulation the 
same augmented van der Waals radii as in the computation of the atomic surface areas. As 
in the previous model the augmented atomic radii, i?^, are set as 

RI = Ri + AR, (26) 

where Ri is the van der Waals radius of the atom and AR = 0.5 A is the offset. The 
augmented radii are used in the same way as in the AGBNPl formulation to define the 
atomic spheres and the associated Gaussian densities [Eqs. ([3])-([6])]. Henceforth in this work 
all of the quantities (atomic volumes, self volumes, etc.) are understood to be computed 
with the augmented atomic radii, unless otherwise specified. In AGBNP2 the form of the 
expression for the inverse Born radii [Eq. 0181) ] is unchanged however the expressions for the 
volume scaling factors Sji and the evaluation of the descreening function Qji are modified as 
follows to introduce the augmented atomic radii and the reduction of the atomic self volumes 
in proportion to the solvent accessible surface RrGclS clS discussed above. 

The pairwise volume scaling factors Sji, that is the volume scaling factor for atom j when 
atom i is removed from the solute, are set as 

Sji = sj + ^ (27) 
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Figure 1: Schematic diagram illustrating the ideas underpinning the model for the solvent excluded 
volume descreening. Circles represent atoms of two idealized solutes placed in proximity of each 
other. The van der Waals description of the molecular volume (panel A) leaves high dielectric 
interstitial spaces that are too small to fit water molecules. The adoption of enlarged van der Waals 
radii (B) removes the interstitial spaces but incorrectly excludes solvent from the surface of solvent 
exposed atoms. The solvent volume subtended by the solvent-exposed surface area is subtracted 
from the enlarged volume of each atom (C) such that larger atomic descreening volumes are assigned 
to buried atoms (circled) than exposed atoms (D), leading to the reduction of interstitial spaces 
while not overly excluding solvent from the surface of solvent exposed atoms. 
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Figure 2: Illustration of the relationship between the van der Waals volume and the solvent 
excluded volume enclosed by the molecular surface. 



where 



(28) 



is the intrinsic scaling factor for atom j (computed with all the atoms present), and the 
quantity 

1 1 1 

(29) 



k k<l 



subtracts from the expression for the self- volume of atom j all those overlap volumes involving 
both atoms i and j. Two differences with respect to the original AGBNPl formulation are 
introduced. The first is that sj [Eq. fl28|) ] is computed from the self-volume after subtracting 
from it the volume of the region subtended by the solvent-exposed surface area in between 
the enlarged and van der Waals atomic spheres of atom j. Referring to Fig. [3] the volume of 
this region is djAj as in Eq. (!28|) with 



3^^ 




(30) 



The other difference concerns the Vj^ term which in the AGBNPl formulation is approxi- 
mated by the 2-body overlap volume Vij] the first term in the right hand side of Eq. (129|) . 
This approximation is found to lack sufficient accuracy for the the present formulation given 
the relative increase in size of all overlap volumes. In the present formulation V^'- is computed 
at all overlap orders according to Eq. (l29il . 

In the AGBNP2 formulation the functional form for the pair descreening function Qji is 
the same as in the original formulation (see Appendix of reference^), however in the new 
formulation this function is evaluated using the van der Waals radius Ri for atom i (the 
atom being "descreened" ) and the augmented radius for atom j (the atom that provides 
the solvent descreening), rather than using the van der Waals radius for both atoms. So if 
the pair descreening function is denoted by Q{r, Ri, R2), where r is the interatomic distance, 
Ri the radius of the atom being descreened and R2 the radius that provides descreening, we 
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dA 




Figure 3: Graphical construction showing the volume subtracted from the atomic self- volume to 
obtain the surface-area corrected atomic self-volume. R is the van der Waals radius of the atom, 
R' = R + Ai? is the enlarged atomic radius. dA is the volume of the region (light gray) subtended 
by the solvent-exposed surface area in between the enlarged and van der Waals atomic spheres. 

set in Eq. ([18]) 

Q,i = Q{rij,Ri,R''^. (31) 

3 Hydrogen Bonding Corrections 

In this section we present the analytical model that implements the short-range hydrogen- 
bonding correction function for AGBNP2. The model is based on a geometrical procedure 
to measure the degree to which a solute atom can interact with hydration sites on the solute 
surface. The procedure is as follows. A sphere of radius Rg representing a water molecule is 
placed in a position that provides near optimal interaction with a hydrogen bonding donor 
or acceptor atom of the solute. The position of this water sphere s is function of the 
positions of two or more parent atoms that compose the functional group including the 
acceptor/donor atom: 

r. = r,({rpj) (32) 

where [vps] represents the positions of the set of parent atoms of the water site s. For 
instance the water site position in correspondence with a polar hydrogen is: 

l^s — I'd + 1 r"HB 

where is the position of the heavy atom donor, vh is the position of the polar hydrogen, 
and d-EB is the distance between the heavy atom donor and the center of the water sphere (see 
Fig. HI). Similar relationships (see Appendix) are employed to place candidate water spheres 
in correspondence with hydrogen bonding acceptor atoms of the solute. These relationships 
are based on the local topology of the hydrogen bonding acceptor group (linear, trigonal, 
and tetrahedral). This scheme places one or two water spheres in correspondence with each 
hydrogen bonding acceptor atom (see Tabled]). 

The magnitude of the hydrogen bonding correction corresponding to each water sphere 
is a function of the predicted water occupancy of the location corresponding to the water 
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Figure 4: Schematic diagram for the placement of the water sphere (w, hght gray) corresponding to 
the hydrogen bonding position relative to the a polar hydrogen (white sphere) of the solute (dark 
gray). The dashed line traces the direction of the hydrogen-parent heavy atom (circled) bond 
along which the water sphere is placed. The magnitude of hydrogen bonding correction grows as a 
function of the volume (light gray) of the water site sphere not occupied by solute atoms. 

sphere. In this work the water occupancy is measured by the fraction Ws of the volume of 
the water site sphere that is accessible to water molecules without causing steric clashes with 
solute atoms (see Fig. (jlj) 

T^free 

= (33) 

where Vg = (4/3)7ri?^ is the volume of the water sphere and 

K''''' = K-E^- + E^-i- E v^^jk (34) 

i i<j i<j<k 

is the "free" volume of water site s, obtained by summing over the two-body, three-body, 
etc. overlap volumes of the water sphere with the solute atoms. Note that the expression 
of the free volume is the same as the expression for the self-volume except that it lacks the 
fractional coefficients 1/2, 1/3, etc. The overlap volumes in Eq. (IMIl are computed using 
radius Rg for the water site sphere (here set to 1.4 A) and augmented radii for the solute 
atoms. Eq. (IMIl is derived similarly to the expression for the self volumes by removing 
overlap volumes from the volume of the water sphere rather than evenly distributing them 
across the atoms participating in the overlap. 

Given the water occupancy Wg of each water sphere, the expression for the hydrogen 
bonding correction for the solute is 

AGhb = '^hsS{ws;Wa, Wb) , (35) 
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Table 1: Optimized surface tension parameters and hydrogen bonding correction parameters 
for the atom types present in protein molecules. 7 is the surface tension parameter, Ny^ is 
the number of water spheres and h is the maximum correction corresponding to each atom 
type [Eq. fl35|) ]. Atom types not listed do not have HB corrections and are assigned 7 = 117 
cal/mol/A^ 

Atom Type 7 (cal/mol/A^) Geometry h (kcal/mol) 



C (aliphatic) 


129 








C (aromatic) 


120 








H on N 




linear 


1 


-0.25 


H on N (Arg) 




linear 


1 


-2.50 


H on 




linear 


1 


-0.40 


H on S 




linear 


1 


-0.50 


(alchol) 


117 


tetrahedral 


2 


-0.40 


(carbonyl) 


117 


trigonal 


2 


-1.25 


(carboxylate) 


40 


trigonal 


2 


-1.80 


N (amine) 


117 


tetrahedral 


1 


-2.00 


N (aromatic) 


117 


trigonal 


1 


-2.00 


S 


117 


tetrahedral 


2 


-0.50 



where hs is the maximum correction energy which depends on the type of solute-water 
hydrogen bond (see Tabled]), and S{w;Wa,Wb) is a polynomial switching function which is 
zero for w < Wa, one for w > Wb, and smoothly (with continuous first derivatives) interpolates 
from to 1 between Wa and Wb (see Fig. [5]). The expression of S{w; Wa, Wb) is: 

{0 w <Wa 

/-(^S) Wa<W<Wb (36) 
1 W > Wb 

where fw{x) is a switching function given by Eq. (Q. In this work we set Wa = 0.15 and 
Wb = 0.5. This scheme establishes (see Fig. [5]) that no correction is applied if more than 75% 
of the water sphere volume is not water accessible, whereas maximum correction is applied 
if 50% or more of the water sphere volume is accessible. 

3.1 Molecular dynamics of mini-proteins 

To be written (usual MD details, time-step, etc.). 

4 Results 

4.1 Accuracy of Born Radii and Surface Areas 

The quality of any implicit solvent model depends primarily on the reliability of the physical 
model on which it is based. The accuracy of the implementation, however, is also a critical 
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Figure 5: The switching function S{w; Wa,Wb) from Eqs. (i36]l and ([9]) with Wa = 0.5 and Wb = 0.9. 

aspect for the success of the model in practice. This is true in particular for models, such 
as AGBNP, based on the Generalized Born formula. It has been pointed out, for instance, 
that approximations in the integration procedure to obtain the Born radii may actually be 
of more significance than the physical approximations on which the GB model is based.— 
It is therefore important to test that the conformational dependent quantities employed by 
AGBNP2 are a good approximation to the geometrical parameters that they are supposed to 
represent. The present AGBNP2 formulation relies mainly on three types of conformational 
dependent quantities: Born radii [Eq. f|T8l) ]. solvent accessible surface areas ffTOj) . and solvent 
accessibilities of hydration sites [Eq. fl36|) ]. In this section we analyze the validity of the 
AGBNP2 analytical estimates for the Born radii and surface areas against accurate numerical 
results for the same quantities. 

We employ the GEPOL program-'^ to compute numerically atomic solvent accessible 
surface areas with a solvent probe diameter of 1 A, the same probe diameter that defines 
the solute-solvent boundary in the AGBNP model. Fig. [6]A shows the comparison between 
the surface area estimates given by the present formulation of AGBNP and the numerical 
surface areas produced by GEPOL for a series of native and modeled protein conformations. 
In Fig.EJB we show the same comparison for the surface areas of the original AGBNPl model. 
These representative results show that the present analytical surface area implementation, 
which as described above employs a weaker switching function for the overlap volumes, 
produces significantly more accurate atomic surface areas than the original model. These 
improvements in the computation of the surface ares, introduced mainly to obtain more 
accurate Born radii through Eq. (1251) . are also expected to yield more reliable cavity hydration 
free energies differences. 

Fig. [7] illustrates on the same set of protein conformations the accuracy of the inverse 
Born radii, B~^, obtained using the AGBNP2 pairwise descreening method using the SEV 
model for the solute volume described above (Eq. (fT8l) ). by comparing them to accurate 
estimates obtained by evaluating the integral in Eq. (IT7|) numerically on a grid. The com- 
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parison is performed for the inverse Born radii because these, being proportional to GB 
self-energies, are a more reliable indicators of accuracy than the Born radii themselves. The 
grid for the numerical integration was prepared as previously reported,^ except that the 
solvent excluding volume (SEV) of the solute was employed here rather than the van der 
Waals volume. The integration grid over the SEV was obtained by taking advantage of 
the particular way that the GEPOL algorithm describes the SEV of the solute; GEPOL 
iteratively places auxiliary spheres of various dimensions in the interstitial spaces between 
solute atoms in such a way that the van der Waals volume of the solute plus the auxiliary 
spheres accurately reproduces the SEV of the solute. Therefore in the present application a 
grid point was considered part of the SEV of the solute if it was contained within any solute 
atom or any one of the auxiliary spheres placed by GEPOL. The default 1.4 A solvent probe 
radius was chosen for the numerical computation of the SEV with the GEPOL program. 

The results of this validation (Fig. [7]) show that the analytical SEV pairwise descreening 
model described above is able to yield Born radii which are not as affected by the spurious 
high dielectric interstitial spaces present in the van der Waals volume description of the 
solute. With the van der Waals volume model (Fig. [7|B) the Born radii of the majority of 
solute atoms start to significantly deviate from the reference values for Born radii larger than 
about 2.5 A {B~^ = 0.4 A~^). Born radii computed with the analytical SEV model instead 
(Fig. [7]A.) track the reference values reasonably well further up to about 4 A {B~^ = 0.25 
A~^). Despite this significant improvement most Born radii are still underestimated by the 
improved model (and, consequently, the inverse Born radii are overestimated - see Fig. [7]), 
particularly those of non-polar atoms near the hydrophobic core of the larger proteins. These 
regions tend to be loosely packed and tend to contain interstitial spaces too large to be 
correctly handled by the present model. Because it mainly involves groups of low polarity 
this limitation has a small effect on the GB solvation energies. It has however a more 
significant effect on the van der Waals solute-solvent interaction energy estimates through 
Eq. fl22|) . which systematically overestimate the magnitude of the interaction for atoms buried 
in hydrophobic protein core. While the present model in general meliorates in all respects 
the original AGBNP model, we are currently exploring ways to address this residual source 
of inaccuracy. 

4.2 Small Molecule Hydration Free Energies 

The validation and parametrization of the hydrogen bonding and cavity correction param- 
eters have been performed based on the agreement between experimental and predicted 
AGBNP2 hydration free energies of a selected set of small molecules, listed in Table [21 con- 
taining the main functional groups present in proteins. This set of molecules includes only 
small and relatively rigid molecules whose hydration free energies can be reliably estimated 
using a single low energy representative conformation— as it was done here. Table [2] lists for 
each molecule the experimental hydration free energy, the AGBNP2 hydration free energy 
computed without HB corrections and the default 7 = 117 cal/mol/A^ surface tension pa- 
rameter, denoted by AGBNP2/SEV, as well as the hydration free energy from the AGBNP2 
model including the HB correction term and the parameters listed in Table [H 

Going down the results in Table [2] we notice a number of issues addressed by the new 
implementation. With the new surface area implementation and without corrections (third 
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Figure 6: Comparisons between numerical and analytical molecular surface areas of the heavy 
atoms of the crystal structures (Ictf and llzl, respectively) of the C-terminal domain of the ribo- 
somal protein L7/L12 (74 aa) and human lysozyme (130 aa), and of four conformations each of 
the trp-cage, cdp-1, and fsd-1 miniproteins extracted from the corresponding explicit solvent MD 
trajectories of the same protein conformations as in Fig. [71 (A) Analytical molecular surface areas 
computed using the present model, and (B), for comparison, analytical surface areas computed 
using the original model from reference.— 
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Table 2: Experimental and predicted hydration free energies of a set of small molecules. 



Molecule 


Exper. 


AGBNP/SEV 


AGBNP2 


n-ethane 


1.83 


0.98 


1.80 


n-propane 


1.9o 


0.92 


1 n'7 

1.97 


n-butane 


2.08 


0.88 


2.14 


n-pentane 


2. 66 


n 7o 
0. ( O 


Z.ZD 


n-nexane 


2.50 


0.70 


2.40 


1 J. 
cyclo-pentane 


1.20 


f\ ^ A 

0.34 


1.63 


cyclo-liexane 


1 OQ 

i.Z6 


U.UO 


i.oU 


benzene 


-0.87 


-1.50 


-1.14 


toluene 


n so 
-U.oy 


-i.oo 


n O/l 
-u.y4 


acetone 


-3.85 


-1.09 


-3.83 


acetoplienone 


-4.0O 


O 7/1 

-Z. 1 4 


e; 07 


ethanol 


-5.01 


A 77 

-4.77 


-5.30 


1 1 

phenol 


-6.62 


A r -1 

-4.51 


-5.38 


ethandiol 


-9.60 


'7 nn 

-7.99 


-9.87 


acetic acid 


-D. ([) 


O 70 

-2. M 


7 O 

- /.05 


propionic acid 


n AO) 

-6.48 


-2.58 


-6.38 


methyl acetate 


-6.61 


-U.IU 


o no 

-o.yz 


ethyl acetate 


-o.iU 


n no 
-U.Uz 


o 

-O.OU 


methyl amine 


-4.00 


o on 

-2.39 


/I 07 

-4.37 


ethyl amine 


-4.oU 


O O/l 

-z.z4 


Q of; 


dimethyl amine 


A C\C\ 

-4.29 


-1.95 


-3.21 


trimethyl amine 


-o.z4 


1 7fi 
-i. ( o 


O QO 

-z.oy 


acetamide 


-9.71 


-6.81 


-10.45 


N-met hylacet amide 


1 n no 

-10.08 


/I Tlx 

-4.75 


7 1 

-7.51 


• 1 * 

pyridine 


-4.70 


-3.62 


-5.30 


r\ J 1 1 • 1 • 

z-methylpyridme 


-4.63 


-2.94 


-4.22 


3- met hy Ipyr idine 


-4. ( I 


-L.oL 


/I 1 Q 

-4. io 


methanethiol 


-1.24 


-0.61 


-1.46 


ethanethiol 


-1.30 


-0.57 


-1.22 


acetate ion 


-79.90 


-77.32 


-87.70 


propionate ion 


-79.10 


-76.29 


-86.29 


methylammonium ion 


-71.30 


-73.21 


-73.54 


ethylammonium ion 


-68.40 


-70.63 


-70.75 


methylguanidinium 


-62.02^ 


-57.30 


-69.81 



In kcal/mol. 

^ Experimental hydration free energy from reference^ except where indicated. 

AGBNP predicted hydration free energies with the default 7 parameter for all atoms types 
(7 = 117 cal/mol/A^) and without HE corrections. 

AGBNP predicted hydration free energies with optimized parameters listed in Table [TJ 
From reference.— 
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Bj'^ (numerical) [A"^] Bf^ (numerical) [A~'] 



Figure 7: Comparisons between numerical and analytical inverse Born radii for the heavy atoms 
of the same protein conformations as in Fig. [6] (A) Analytical Born radii computed using the 
present SEV model. (B) Analytical Born radii computed using the van der Waals volume model 
(reference^). 

column in Table [2]) the hydration free energies of the normal alkanes are too small compared 
to experiments, furthermore, in contrast with the experimental behavior, the predicted hy- 
dration free energies incorrectly become more favorable with increasing chain length. A 
similar behavior is observed for the aromatic hydrocarbons. Clearly this is due to the rate of 
increase of the positive cavity term with increasing alkane size which is insufficient to offset 
the solute-solvent van der Waals interaction energy term, which becomes more negative with 
increasing solute size. We have chosen to address this shortcoming by increasing by 10.2% 
and 2.5% respectively the effective surface tensions for aliphatic and aromatic carbon atoms 
rather than decreasing the corresponding a parameters of the van der Waals term since the 
latter had been previously validated against explicit solvent simulations. We have chosen to 
limit the increases of the surface tension parameters to aliphatic and aromatic carbon atoms 
since the results for polar functional groups did not indicate that this change was neces- 
sary for the remaining atom types. With this new parametrization we achieve (compare the 
second and fourth column in Table [2]) excellent agreement between the experimental and 
predicted hydration free energies of the alkanes and aromatic compounds. Note that the 
AGBNP2 model, regardless of the parametrization, correctly predicts the more favorable 
hydration free energies of the cyclic alkanes relative to their linear analogous. AGBNP2, 
thanks to its unique decomposition of the non-polar solvation free energy into unfavorable 
cavity term and an opposing favorable term, is, to our knowledge, the only analytic implicit 
solvent implementation capable of describing correctly this feature of the thermodynamics 
of hydration of hydrophobic solutes. 

The AGBNP2 model without corrections markedly underpredicts the magnitudes of the 
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experimental hydration free energies of the compounds containing carbonyl groups (ketones, 
organic acids, and esters). The hydration free energies of alcohols are also underpredicted but 
by smaller amounts. Much better agreement with the experimental hydration free energies of 
these oxygen-containing compounds (see Table [2]) is achieved after applying hydrogen bond- 
ing corrections with h = —1.25 kcal/mol for the carbonyl oxygen, and h = —0.4 kcal/mol for 
both the hydrogen and oxygen atoms of the hydroxy group (Tabled]). Note that the same pa- 
rameter employed individually for carbonyl and hydroxy groups in ketones and alcohols are 
applied to the more complex carboxylic groups of acids and esters as well as amides and car- 
boxylate ions. The thiol group of organic sulfides required similar corrections as the hydroxy 
group (Tables [H and [21). The AGBNP2 model without corrections also markedly underpre- 
dicted the magnitude of the experimental hydration free energies of amines and amides, and, 
to a smaller extent, of compounds with nitrogen containing etherocyclic aromatic rings. The 
addition of HB corrections of —0.25 kcal/mol for amine hydrogens and h = —2.0 kcal/mol for 
both amine and aromatic nitrogen atoms yields improvement agreement (Table [2]), although 
the effect of N-methylation is still over-emphasized by both AGBNP2 parametrization. 

4.3 Mini-protein results 

As described in Section 13.11 we have performed restricted MD simulations of a series of 
so-called mini-proteins (trp-cage, cdp-1, and fsd-1) to study the extent of the agreement 
between the conformational ensembles generated with the original AGBNP implementation 
(AGBNPl) and the present implementation (AGBNP2) with respect to explicit solvent- 
generated ensembles. Earlier studies suggested that the AGBNP/OPLS-AA model correctly 
reproduced for the most part backbone secondary structure features of protein and peptides. 
The tests in the present study are therefore focused on sidechain conformations. The back- 
bone atoms were harmonically restricted to remain within approximately 3 A C-a RMSD 
of the corresponding NMR experimental models. We structurally analyzed the ensembles in 
terms of the extent of occurrence of intramolecular interactions. 

As shown in Table [3] we measured a significantly higher average number of intramolecular 
hydrogen bonds and ion pairing in the AGBNPl ensembles relative to the explicit solvent 
ensembles for all mini-proteins studied. The largest deviations are observed for cdp-1 and 
fsd-1, two mini-proteins particularly rich in charged sidechains, with on average nearly twice 
as many intramolecular hydrogen bonds compared to explicit solvent. Many of the excess 
intramolecular hydrogen bonds with AGBNPl involve interactions between polar groups 
(polar sidechains or the peptide backbone) and the sidechains of charged residues. For 
example for cdp-1 we observe (see Table [3]) approximately 8 hydrogen bonds between polar 
and charged groups on average compared to nearly none with explicit solvation. 

Despite the introduction of empirical surface area corrections to penalize ion pairs,— 
AGBNPl over-predicts ion-pair formation. We found that ion pairing involving arginine 
were particularly over-stabilized by AGBNPl as we observed stable ion pairing between 
arginine and either glutamate or aspartate residue during almost the entire duration of the 
simulation in virtually all cases in which this was topologically feasible given the imposed 
backbone restrains. In contrast, with explicit solvation some of the same ion pairs were 
seen to form and break numerous times indicating a balanced equilibrium between contact 
and solvent-separated conformations. This balance is not reproduced with implicit solvation 
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Table 3: Average number of some types of intramolecular electrostatic interactions in the 
explicit solvent conformational ensembles, and the AGBNPl and AGBNP2 ensembles of the 
trp-cage, cdp-1, and fsd-1 mini-proteins. 

3lini-protein Explicit AGBNPl AGBNP2 
Intramolecular Hydrogen Bonds 
trp-cage 13.5 18.3 15.3 

cdp-1 12.6 24.5 15.4 

fsd-1 14.1 24.6 14.3 

Polar-Polar Hydrogen Bonds 
trp-cage 12.9 17.1 13.8 

cdp-1 12.5 16.4 14.1 

fsd-1 12.0 14.2 12.9 



Polar-Charged Hydrogen Bonds 



trp-cage 


0.6 


1.2 


1.4 


cdp-1 


0.1 


8.1 


1.3 


fsd-1 


2.1 


9.6 


1.4 




Ion 


Pairs 




trp-cage 


0.3 


1.0 


1.0 


cdp-1 


2.5 


2.9 


2.7 


fsd-1 


1.4 


4.6 


4.0 



which instead strongly favors ion pairing. In any case, the relative stability of ion pairs 
appeared to depend in subtle ways on the protein environment as for example the two ion 
pairs between arginine and glutamate of cdp-1 were found to be stable with either explicit 
solvation or AGBNPl implicit solvation whereas other Arg-Glu ion pairs in trp-cage and 
fsd-1 were found to be stable only with implicit solvation. 

This analysis generally confirms quantitatively a series of past observations made in our 
laboratory indicating that the original AGBNP implementation tends to be biased towards 
conformations with excessive intramolecular electrostatic interactions, at the expense of more 
hydrated conformations in which polar groups are oriented so as to interact with the water 
solvent. During the process of development of the modifications to address these problems, 
we found useful to rescore with varying AGBNP2 formulations and parametrization the mini- 
protein conformational ensembles obtained with AGBNPl and explicit solvation, bypassing 
the time consuming MD runs to obtain ensembles with each new parametrization. An ex- 
ample of this analysis is shown in Fig. [HI which compares the probability distributions of the 
AGBNPl and AGBNP2 effective potential energies over the conformational ensembles gener- 
ated with AGBNPl implicit solvation and with explicit solvation. These results clearly show 
that the AGBNPl/OPLS-AA effective potential disfavor conformations from the explicit 
solvent ensemble relatively to those generated with implicit solvation. The AGBNPl energy 
scores (first row of Fig. [8]) of the explicit solvent ensembles of all mini-proteins are shifted 
towards higher energies than those of the AGBNPl ensemble, indicating that conformations 
present in the explicit solvent ensemble would be rarely visited when performing confor- 
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mational sampling with the AGBNPl/OPLS-AA potential. AGBNPl/OPLS-AA assigns a 
substantial energetic penalty (see Fig. [8]A-C) to the explicit solvent ensemble relative to the 
AGBNPl ensemble (on average 3.3, 4.4, and 5.7 kcal/mol per residue for, respectively, the 
trp-cage, cdp-1, and fsd-1 mini-proteins). This energetic penalty, being significantly larger 
than thermal energy, rules out the possibility that conformational entropy effects could off- 
set it to such an extent so as to equalize the relative free energies of the two ensembles. 
Detailed analysis of the energy scores shows that, as expected, the AGBNPl implicit solvent 
ensemble is mainly favored by more favorable electrostatic Coulomb interaction energies due 
to its greater number of intramolecular electrostatic contacts relative to the explicit solvent 
ensemble (see above). Conversely, the AGBNPl solvation model does not assign sufficiently 
favorable hydration free energy to the more solvent-exposed conformations obtained in ex- 
plicit solvent so as to make them competitive with the more compact conformations of the 
AGBNPl ensemble. 

The energy distributions of the AGBNP2/0PLS-AA effective energy function for the 
same ensembles (the AGBNPl ensemble and the explicit solvent ensemble) showed marked 
improvements (see Fig. [HI second and third row). The second row of graphs in Fig. [8] (panels 
D-F) correspond to the AGBNP2 model without hydrogen bonding to solvent corrections 
(referred to as AGBNP2-SEV). We see that the introduction of the SEV model for the solute 
volume already drastically lowers the energy scores of the explicit solvent conformations 
relative to the implicit solvent conformations. The AGBNP2-SEV energy gaps between the 
explicit and implicit solvation ensembles are reduced approximately by a factor of two for all 
three mini-proteins, moreover the energy distributions for the two ensembles now partially 
overlap. The AGBNP2/0PLS-AA energy scores with the hydrogen bonding corrections and 
the parameters in Table [1] (Figs. [HlG-I) favor the explicit solvent ensemble relative to the 
AGBNPl ensemble even further. The average energy gaps between the two ensembles are 
reduced to 1.3, 1.0, and 1.2 kcal/mol per residue for, respectively, the trp-cage, cdp-1, and fsd- 
1 mini-proteins, compared to the average energy gaps of 3.3, 4.4, and 5.7 kcal/mol with the 
AGBNPl/OPLS-AA effective potential (a 3 to 5-fold improvement) and the overlaps between 
the energy distributions of the implicit and explicit solvation ensembles is further increased 
relative to AGBNP2-SEV (compare Figs.[Hn-I with Figs.[HP-F). The energy bias per residue 
of AGBNP2 relative to the explicit solvent ensemble is comparable to thermal energy and 
smaller than the spread of the energy distributions indicating that AGBNP2 is capable of 
generating explicit solvent-like conformations much more frequently than AGBNPl. 

The energy scoring experiments on the explicit solvent and AGBNPl ensembles described 
above were very useful for tuning the formulation of the AGBNP2 model without requiring 
running lengthy MD simulations. They do not, however, guarantee that the conformational 
ensembles generated with the AGBNP2 solvation model will more closely match the explicit 
solvent ensembles than those generated with AGBNPl. This is because the new solvation 
model could introduce new energy minima not encountered with AGBNPl or explicit solva- 
tion that would be visited only by performing conformational sampling with AGBNP2. To 
validate the model in this respect we obtained MD trajectories with the AGBNP2 implicit 
solvent model and compared the corresponding probability distributions of the effective en- 
ergy with those of the explicit solvent ensembles similarly as above. The results for the 
three mini-proteins, shown in Fig. [H indicate that the AGBNP2-generated ensembles dis- 
play significantly smaller bias (mean energy gaps per residue of 2.0, 2.1, and 2.5 kcal/mol 
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Figure 8: Probability distributions of the AGBNPl/OPLS-AA (first row, panels A, B, and C), 
AGBNP2-SEV/0PLS-AA (second row, panels D, E, and F), and AGBNP2/0PLS-AA effective 
potential energies for the conformational ensembles for the trp-cage (first column, panels A, D, 
G), cdp-1 (second column, panels B, E, H), and fsd-1 (third column, panels C, F, I) mini-proteins 
obtained using the AGBNPl solvation model (full line) and explicit solvation (dashed line). The 
AGBNP2-SEV model is the AGBNP2 model with the SEV representation of the solute volume and 
without the hydrogen bonding corrections [hs = in Eq. (f35l) for all atom types]. The distributions 
are shown as a function of the energy gap per residue (Au) relative to the mean effective potential 
energy of the implicit solvent ensemble. 
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Figure 9: Probability distributions of the AGBNP2/0PLS-AA effective potential energy for the 
conformational ensembles of the trp-cage (A), cdp-1 (B), and fsd-1 (C) mini-proteins obtained 
using the present AGBNP2 solvation model (full line) and explicit solvation (dashed line). The 
distributions are shown versus the energy gap per residue (Au) relative to the mean effective 
potential energy of the implicit solvent ensemble. 

for, respectively, the trp-cage, cdp-1, and fsd-1 mini-proteins) than AGBNPl (Fig. [HI panels 
A-C) which instead yielded energy gaps of 3.3, 4.4, and 5.7 kcal/mol per residue respectively. 
This observation confirms that conformational sampling with AGBNP2 produces conforma- 
tions that more closely match on average the explicit solvent ensemble without producing 
unphysical minima that deviate significantly from it. 

We have structurally analyzed the conformational ensembles obtained with the AGBNPl 
and AGBNP2 models to establish the degree of improvement achieved with the new model 
in terms of occurrence of intramolecular interactions. The salient results of this analysis are 
shown in Table [3l This table reports for each mini-protein the average number of intramolec- 
ular hydrogen bonds and ion pairs. The number of hydrogen bonds is further decomposed 
into those involving only polar groups (including the backbone) and those involving a polar 
groups and the sidechain of a charged residue (arginine, lysine, aspartate, and glutamate). 
As noted above it is apparent from this data that the AGBNPl model produces conforma- 
tions with too many hydrogen bonds and ion pairs. The majority of the excess hydrogen 
bonds with AGBNPl involve aminoacid sidechains. Similarly, too many ion pairs are ob- 
served in the AGBNPl ensemble particularly for the fsd-1 mini-protein (4.6 ion pairs on 
average with AGBNPl compared to only 1.4 in explicit solvent). The AGBNP2 ensembles, 
in comparison, yields considerably fewer intramolecular hydrogen bonds. For instance, the 
average number of hydrogen bonds for cdp-1 is reduced from 24.5 with AGBNPl to 15.4 
with AGBNP2. With AGBNP2 the number of polar-polar hydrogen bonds is generally in 
good agreement with explicit solvation. However the greatest improvement is observed with 
polar-charged interactions. For example the number of polar-charged hydrogen bonds of 
fsd-1 is reduced by almost ten fold in going from AGBNPl to AGBNP2 to reach good agree- 
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ment with explicit solvation. Importantly, a significant fraction of the excess polar-charged 
interactions observed with AGBNPl corrected by AGBNP2 are interactions between the 
peptide backbone and charged sidechains that would otherwise interfere with the formation 
of secondary structures. 

With AGBNP2 we observe small but promising improvements in terms of ion pair forma- 
tion. The average number of ion pairs of cdp-1 consistently agrees between all three solvation 
models, and the only possible ion pair in trp-cage is more stable in both implicit solvent for- 
mulations than in explicit solvent (it occurs in virtually all implicit solvent conformations 
compared to only 30% of the conformations in explicit solvent). However the average number 
of ion pairs for fsd-1 is reduced from 4.6 with AGBNPl to 4.0 with AGBNP2. We observe 
good agreement between the number of ion pairs involving lysine with either AGBNPl or 
AGBNP2 and explicit solvation. However ion pairs involving arginine are generally more 
stable with implicit solvation than with explicit solvation. The agreement in the number of 
ion pairs with cdp-1 is due to the fact that for this mini-protein the two possible ion pairs 
involving arginine result stable with explicit solvation as well as with implicit solvation. For 
the other two mini-proteins however, ion pairs involving arginine that are marginally sta- 
ble with explicit solvation are found to be significantly more stable with implicit solvation, 
although less so with AGBNP2 solvation. 

5 Discussion 

The non-linear and asymmetric dielectric response of water stems primarily by the finite 
extent and internal structure of water molecules.- As further discussed below, the model- 
ing of effects due to water granularity is important for the proper description of molecular 
association equilibria. Integral equation methods^ provide an accurate implicit solvation 
description from first principles, however, despite recent progress,— they are not yet appli- 
cable to molecular dynamics of biomolecules. The primary aim of the present study has been 
to formulate an analytical and computational efficient implicit solvent model incorporating 
solvation effects beyond those inherent in standard continuum dielectric models and, by so 
doing, achieving an improved description of solute conformational equilibria. 

In this work we have developed the AGBNP2 implicit solvent model which is based on 
an empirical (but physically-motivated) first solvation shell correction function parametrized 
against experimental hydration free energies of small molecules and the results of explicit 
solvent molecular dynamics simulations of a series of mini-proteins. The correction function 
favors conformations of the solute in which polar groups are oriented so as to form hydrogen 
bonds with the surrounding water solvent thereby achieving a more balanced equilibrium 
with respect to the competing intramolecular hydrogen bond interactions. A key ingredient 
of the model is an analytical prescription to identify and measure the volume of hydration 
sites on the solute surface. Hydration sites that are deemed too small to contain a water 
molecule do not contribute to the solute hydration free energy. Conversely, hydration sites 
of sufficient size form favorable interactions with nearby polar groups. This model thus 
incorporates the effects of both water granularity and of non-linear first shell hydration 
effects. 

The GB and non-polar models in the AGBNP2 implicit solvent model provide the linear 
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continuum dielectric model basis of the model as well as a description of non-electrostatic hy- 
dration effects.— In this work the GB and solute-solvent dispersion interaction energy models 
are further enhanced by replacing the original van der Waals solute volume model with a 
more realistic solvent excluding volume (SEV) model. The new volume description improves 
the quality of the Born radii of buried atoms and atoms participating in intramolecular in- 
teractions which would otherwise be underestimated due to high dielectric interstitial spaces 
present with the van der Waals volume description.— GB models with these characteristics 
have been previously proposed. The GBMV series of models^i^i^ represent the SEV on a 
grid which, although accurate, is computationally costly and lacks frame of reference invari- 
ance. The pairwise descreening-based GB^^*-^ model-l^I introduced an empirical rectifying 
function to increase the Born radii of buried atoms in an averaged, geometry-independent 
manner. The GBn model— introduced the neck region between pairs of atoms as additional 
source of descreening, dampened by empirical parameters to account in an average way for 
overlaps between neck regions and between solute atoms and neck regions. The approach 
proposed here to represent the SEV consists of computing the atomic self- volumes [Eq. (IT^ ]. 
used in the pairwise descreening computation, using enlarged atomic radii so as to cover the 
interatomic interstitial spaces. The self-volume of each atom is then reduced proportionally 
to its solvent accessible surface area [see Eq. fl28l) ] to subtract the volume in van der Waals 
contact with the solvent. We show (Fig. [7l) that this model reproduces well Born radii com- 
puted from an accurate numerical representation of the SEV, noting that improvements for 
the Born radii of atoms in loosely packed hydrophobic interior, while significant, are still 
not optimal. Although approximate, this representation of the SEV maintains the simplicity 
and computational efficiency of pairwise descreening schemes, while accounting for atomic 
overlaps in a consistent and parameter-free manner. 

The new AGBNP2 model has been formulated to be employed in molecular dynamics 
conformational sampling applications, which require potential models of low computational 
complexity and favorable scaling characteristics, and with analytical gradients. These re- 
quirements have posed stringent constraints on the design of the model and the choice of 
the implementation algorithms. In the formulation of AGBNP2 we have reused as much as 
possible well-established and efficient algorithms developed earlier for the AGBNPl model. 
For example the key ingredient of the hydrogen bonding correction function is the free vol- 
ume of an hydration site, which is computed using a methodology developed for AGBNPl 
to compute atomic self-volumes. Similarly the SEV-based pairwise descreening procedure 
employs atomic surface areas (see Eq. [28]), computed as previously described."^ AGBNP2 
suffers additional computational cost associated with the SEV-based pairwise descreening 
procedure and the hydrogen bonding correction function. This handicap however is offset by 
having only one solute volume model in AGBNP2 rather than two in AGBNP2. AGBNPl 
requires two separate invocations of the volume overlaps machinery [Eq. [2] for each of the 
two volume models it employs, corresponding to the van der Waals atomic radii for the pair- 
wise descreening calculation and enlarged radii for the surface area calculation.^^ AGBNP2 
instead employs a single volume model for both the pairwise descreening and surface area 
calculations. A direct CPU timing comparison between the two models can not be re- 
ported at this time because the preliminary implementation of the AGBNP2 computer code 
used in this work lacks key data caching optimizations similar to those already employed in 
AGBNPl. Given the computational advantages of the new model discussed above we expect 
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to eventually obtain similar or better performance than with AGBNPl. 

The AGBNP2 model has been parametrized against experimental hydration free energies 
of a series of small molecules and with respect to the ability of reproducing energetic and 
structural signatures of the conformational ensemble of three mini-proteins generated with 
explicit solvation. These data sources are chosen so as to ensure that the resulting model 
would be applicable to both hydration free energy estimation and conformation equilibria, 
which are fundamental characteristics for models aimed at protein-ligand binding affinity 
estimation. On the other hand, experimental hydration free energies and explicit solvent 
conformational ensembles are to some extent incompatible with one another given the limi- 
tations of even the best fixed charge force fields and explicit solvation models to reproduce 
experimental hydration free energies of small molecules with high accuracy.— i^*^ Mind- 
ful of this issue we did not seek a perfect correspondence with the experimental hydration 
free energy values. We first obtained parameters by fitting against the small molecule ex- 
perimental hydration free energies, then adjusted the parameters to improve the agreement 
with the explicit solvent data making sure that the predicted small molecule hydration free 
energies remained within a reasonable range relative to the experimental values. In practice 
this procedure yielded predicted hydration free energies in good agreement with the exper- 
imental values with the exception of the carboxylate and guanidinium ions (see Table [2]) 
for which AGBNP2 predicts more favorable hydration free energies than the experiments; a 
consequence of the large hydrogen bonding corrections necessary to reduce the occurrence 
of intramolecular electrostatic interactions in the investigated proteins. As discussed further 
below, limitations in the description of hydration sites adopted for carboxylate and guani- 
dinium ions may be partly the cause of the observed inconsistencies for these functional 
groups. 

Parametrization and validation of the model focused in part on comparing the effective 
potential energy distributions of implicit solvent conformational ensembles with those of ex- 
plicit solvent ensembles. We observed that the AGBNPl solvation model energetically ranked 
explicit solvent conformations significantly less favorably than implicit solvent conforma- 
tions. The AGBNP2 model is characterized by smaller energetic bias relative to the explicit 
solvent ensembles, indicating that conformational sampling with the AGBNP2/0PLS-AA 
energy function produces conformations that more loosely match those obtained with ex- 
plicit solvation. This result is a direct consequence of employing the more realistic solvent 
excluding volume description of the solute which yields larger Born radii for buried groups, as 
well as the hydrogen bonding to solvent corrections, which favor solvent exposed conforma- 
tions of polar groups. Furthermore, comparison of the energy distributions of the AGBNP2 
and explicit solvent ensembles for the three mini-proteins (Fig. ^ shows, in contrast to the 
AGBNPl results, that the AGBNP2 bias for the two more charge-rich mini-proteins (cdp-1 
and fsd-1) is similar to that of the least charged one (trp-cage). This suggests that the 
residual energetic bias of the AGBNP2 model is probably related to the non-polar model 
rather than the electrostatic model. Future studies will address this particular aspect of the 
model. 

The energy scoring studies conducted in this work indicate that AGBNP2 is a significant 
improvement over AGBNPl. They also show, however, that the new model falls short of 
consistently scoring explicit solvent conformations similarly to implicit solvent conforma- 
tions. Although an optimal match between energy distributions is a necessary condition for 
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complete agreement between implicit and explicit solvation results, it is unrealistic to expect 
to reach this ultimate goal at the present level of model simplification. Increasing the mag- 
nitude of the hydrogen bonding corrections can improve the agreement between the explicit 
and implicit solvation energy distributions, albeit at the expense of the quality of the pre- 
dicted small molecule hydration free energies. It seems likely that the no parametrization of 
the current model would yield both good relative conformational free energies and hydration 
free energies. Future work will pursue the modeling of additional physical and geometrical 
features necessary to improve the agreement between implicit and explicit solvation energy 
distributions. The energy gap between the implicit solvent and explicit solvent energy dis- 
tributions used here is, we believe, a meaningful measure of model quality and could serve 
as a useful general validation tool to compare the accuracy of implicit solvent models. 

The excessive number of intramolecular electrostatic interactions involving charged groups 
has been one of the most noticeable shortcoming of GB-based implicit solvent models.— 
To correct this tendency we have in the past adopted in the AGBNPl model ad-hoc strate- 
gies aimed at either destabilizing electrostatic intramolecular interactions^ or, alternatively, 
stabilizing the competing solvent-separated conformations.— This work follows the latter 
approach using a more robust and physically-motivated framework based on locating and 
scoring hydration sites on the solute surface as well as adopting a more realistic volume 
model. Structural characterization of the conformational ensembles has shown that AGBNP2 
produces significantly fewer intramolecular interactions than AGBNPl reaching good agree- 
ment with explicit solvent results. The reduction of intramolecular interactions has been the 
greatest for interactions between polar and charged groups. We believe that the excessive 
tendency toward the formation of intramolecular interactions to be the root cause of the 
high melting temperatures of structured peptides^ predicted with AGBNPl. Given the 
reduction of intramolecular interactions achieved with AGBNP2, we expect the new model 
to yield more reasonable peptide melting temperatures; a result which we hope to report in 
future publications. 

Less visible improvements have been obtained for ion pairs involving arginine sidechains 
which remain more stable with implicit solvation than explicit solvation. However, signif- 
icantly, with AGBNP2 we observed a more dynamical equilibrium between ion-paired and 
solvent-separated conformations of arginine sidechains that was not observed with AGBNPl. 
This result is promising because it indicates that the AGBNP2 solvation model, although 
still favoring ion paired conformations, produces a more balanced equilibrium, which is in- 
stead almost completely shifted towards contact conformations with AGBNPl. Nevertheless 
it is apparent that the AGBNP treatment of the guanidinium group of arginine is not as 
good as for other groups. This limitation appears to be shared with other functional groups 
containing sp2-hybridized nitrogen atoms as evidenced for example by the relatively lower 
quality of the hydration free energy predictions for amides and nitrogen-containing aro- 
matic compounds (Table [2]). Similar implicit solvent overstabilization solvation of arginine- 
containing ion pairs has been observed Yu et al.~ in their comparison of Surface Generalized 
Born (SGB) and SPG explicit solvation with the OPLS-AA force field. Despite quantita- 
tive differences, the explicit solvent studies (with the TIP3P water model) of Masunov and 
Lazaridis^ and Hassan,— using the GHARMM force field, and that of Mandell at al.,— 
using the OPLS-AA force field, have confirmed that arginine forms the strongest ion pairing 
interactions, especially in the bidentate coplanar conformation. These observations are con- 
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sistent with the present exphcit solvent results using OPLS-AA and the SPC water model, 
where we find that most of the ion pairs of the mini-proteins were found to involve arginine 
sidechains. In contrast to our present implicit solvent results, however, the work of Masunov 
and Lazaridis"^^ indicated that the GB-based implicit solvent model that they analyzed-^^ 
produced potentials of mean force for arginine-containing ion pairs in general agreement with 
explicit solvation. 

To rationalize the present implicit solvent results concerning ion pair formation, it has 
been instructive to analyze the potentials of mean force (PMF) of ion pair association with 
the AGBNP model. As an example. Fig. shows the PMF for the association of propyl- 
guanidinium (arginine sidechain analog) and ethyl-acetate (aspartate and glutamate analog) 
in a bidentate coplanar conformation (similar to the arrangement used previously— »^>^>^) 
for various AGBNP implementations. The corresponding explicit solvent PMF obtained by 
Mandell et al.— is also shown in Fig. [TO]for comparison. The original AGBNPl parametriza- 
tion^' clearly leads to an overly stable salt bridge with the contact conformation scored at 
approximately —19 kcal/mol relative to the separated conformation, compared with —8.5 
kcal/mol with explicit solvation. The AGBNPl parametrization analyzed here, which in- 
cludes an empirical surface area correction to reduce the occurrence of ion pairs,— yields 
a contact free energy (—11 kcal/mol) in much better agreement with explicit solvation, al- 
though the shape of the PMF is not properly reproduced. The present AGBNP2 model 
without hydrogen bonding corrections (labeled "AGBNP2-SEV" in Fig. [10]) yields a PMF 
intermediate between the original and corrected AGBNPl parametrization. The AGBNP2 
model with hydrogen bonding corrections yields the PMF with the closest similarity with the 
one obtained in explicit solvent (Fig. [TUB). Not only the contact free energy (—6.5 kcal/mol) 
is in good agreement with the explicit solvent result, but, importantly, it also reproduces the 
solvation barrier of the PMF at 5 A separation, corresponding to the distance below which 
there is insufficient space for a water layer between the ions. 

It is in this range of distances that the greatest discrepancies are observed between PMF's 
with explicit solvation and some GB-based implicit solvation models^-'-^- that do not model 
effects due to the finite size of water molecules. Both the hydrogen bonding correction and 
the SEV volume description employed in AGBNP2, which are designed to take into account 
the granularity of the water solvent - the hydrogen bonding correction through the minimum 
free volume of water sites [Eq. ([35!) ] and the SEV model through the water radius offset [Eq. 
( !26|) ] - make it possible to reproduce the solvation barrier typical of molecular association 
processes in water. 

The good correspondence between the AGBNP2 and explicit solvent PMF's for propyl- 
guanidinium and ethyl-acetate (Fig. [TUl) stand in contrast with the residual AGBNP2 over- 
prediction of arginine salt bridges compared to explicit solvation (Table [3]). We observed 
that in the majority of arginine salt bridges occurring with AGBNP2, the guanidinium and 
carboxylate groups interact at an angle rather than in the coplanar configuration discussed 
above. We have confirmed that the PMF of ion pair formation for an angled conformation 
(not shown) indeed shows a significantly more attractive contact free energy than the copla- 
nar one. This result indicates that the in-plane placement of the hydration sites for the car- 
boxylate groups (see Appendix) does not sufficiently penalize angled ion pair arrangements. 
This observation is consistent with the need of introducing of an isotropic surface-area based 
hydration correction for carboxylate groups (the reduced 7 parameter for the carboxylate 
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Figure 10: Potential of mean force of ion pair formation between propyl-guanidinium and ethyl- 
acetate in the coplanar orientation with AGBNP implicit solvation (A) and explicit solvation 
(B; referenc&ii^) . In (A) "AGBNPl (orig.)" refers to the original AGBNPl parametrization,— 
"AGBNPl" refers to the AGBNPl model used in this work which includes a surface tension pa- 
rameter correction for the carboxylate group aimed at reducing the occurrence of ion pairs,— 
"AGBNP2" refers to the current model and "AGBNP2-SEV" the current model without hydrogen 
bonding corrections. The potentials of mean force are plotted with respect to the distance between 
the atoms of the protein sidechain analogs equivalent to the of arginine and the C7 of aspartate. 
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oxygen atoms), which showed some advantage in terms of reducing the occurrence of salt 
bridges. Future work will focus on developing a more general hydration shell description for 
carbonyl groups and related planar polar groups to address this issue. 

6 Conclusions 

We have presented the AGBNP2 implicit solvent model, an evolution of the AGBNPl model 
we have previously reported, with the aim of incorporating modeling of hydration effects 
beyond the continuum dielectric representation. To this end a new hydration free energy 
component based on a procedure to locate and score hydration sites on the solute surface 
is used to model first solvation shell effects, such as hydrogen bonding, which are poorly 
described by continuum dielectric models. This new component is added to the Generalized 
Born and non-polar AGBNP models which have been improved with respect to the descrip- 
tion of the solute volume description. We have introduced an analytical Solvent Excluding 
Volume (SEV) model which reduces the effect of artefactual high-dielectric interstitial spaces 
present in conventional van der Waals representations of the solute volume. The new model 
is parametrized and tested with respect to experimental hydration free energies and the 
results of explicit solvent simulations. The modeling of the granularity of water is one of 
the main principles employed in the design of the empirical first shell solvation function and 
the SEV model, by requiring that hydration sites have a minimum available volume based 
on the size of a water molecule. We show that the new volumetric model produces Born 
radii and surface areas in good agreement with accurate numerical evaluations. The results 
of Molecular Dynamics simulations of a series of mini-proteins show that the new model 
produces conformational ensembles in much better agreement with reference explicit solvent 
ensembles than the AGBNPl model both with respect to structural and energetics measures. 

Future development work will focus on improving the modeling of some functional groups, 
particularly ionic groups involving sp2 nitrogen, which we think are at the basis of the residual 
excess occurrence of salt bridges, and on the optimization of the AGBNP2 computer code 
implementation. We also plan to test the model on wide variety of benchmarks including 
peptide folding. 

A Hydration Site Locations 

Fig. [11] shows the location of the hydration sites for the functional groups listed in Table [H 
Each hydration site is represented by a sphere of 1.4 A radius. The distance cinB between 
the donor or acceptor heavy atom and the center of the hydration site sphere is set to 2.5 A. 

There is a single linear geometry for HB donor groups. The corresponding hydration site 
is placed at a distance cinB from the heavy atom donor along the heavy atom-hydrogen bond. 

Acceptor trigonal geometries have one or two hydration sites depending on whether the 
acceptor atom is bonded to, respectively, two or one other atom. In the former case the 
water site is placed along the direction given by the sum of the unit vectors corresponding 
to the sum of the NRi and NR2 bonds (following the atom names in Fig. [TT]) . In the latter 
case the Wi site (see Fig. [TTj) is placed in the RiCO plane forming an angle of 120° with the 
CO bond. The W2 site is placed similarly. 
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Acceptor tetrahedral geometries have one or two hydration sites depending on whether 
the acceptor atom is bonded, respectively, to three or two other atoms. In the former case 
the water site is placed along the direction given by the sum of the unit vectors corresponding 
to the sum NRi, NR2, and NR3 bonds. In the latter case the positions of the Wi and W2 
sites are given by 

wi = O + duB (cos ^ui + sin 9u2) 
wi = O + duB (cos 9ui — sin 9u2) 

where O is the position of the acceptor atoms, 6 = 104.4°, and ui and U2 are, respectively, 
the unit vectors corresponding to the ORi and OR2 bonds. 



B Gradients of GB and van der Waals Energies 

The component of the gradient of the AGBNP2 van der Waals energy at constant self- 
volumes is the same as in the AGBNPl model [see Appendix C of reference^]. In AGBNP2 
the expression for the component of the gradient corresponding to variations in the atomic 
scaling factors, Sij, includes pair corrections at all overlap levels because of the presence of 
multi-body volumes in V-j. In addition, a new component corresponding to the change in 
surface areas appears: 

= -—T-^Qki (37) 

Eq. ( I37I) leads to the same expression of the derivative component as in the AGBNPl model 
[Eq. (72) in reference^] (except for the extra elements in the 2-body terms due to the 
inclusion of the 1/2 V^j correction term). Eq. fl38|) corresponds to the component of the 
derivative due to variations in Vjf., the volume to be added to the self- volumes of j and k 
to obtain the sjk and Skj scaling factors. In the AGBNPl model this component included 
only 2-body overlap volumes; in AGBNP2 this term instead includes overlap volumes at all 
orders. Finally, Eq. ( 139|) . where Ak is the surface area of atom k, leads to the component 
of the derivatives of the GB and vdW terms due to variations of the exposed surface area. 
The latter two terms are new for AGBNP2. 



B.l Component of derivative from Eq. (1381) 

From Eq. (63) in reference^ and Eq. fl38l) we have 



4- =EW^^.^, (40) 

\ / Q2 jk '^'■i 
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Tetrahedral(3) Tetrahedral(2) 



Figure 11: Diagram illustrating the hydration site locations for each of the functional group 
geometries used in this work. Linear: hydrogen bond donor; trigonal(l) and trigonal(2): 
trigonal planar geometries with, respectively, one and two covalent bonds on the acceptor 
atom; tetrahedral(2) and tetrahedral(3): tetrahedral geometries with, respectively, two and 
three covalent bonds on the acceptor atom. Representative molecular structures are shown 
for each geometry. 
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where Wkj has the same expression as in Eq. (69) in reference.— In working with Eq. (HOl) 
it is important to note that, whereas V^j is symmetric with respect to swapping the j and k 
indexes, Wkj and Wjk are different from each other. Substituting Eq. (l29ll into Eq. ( I40l) and 
expanding over symmetric terms we obtain 



(41) 



Eq. (HTj) is simphfied by noting that 



^ = + fo^ + , , , (42) 

OTj OTj OTj 

Eq. fH2l) is inserted in Eq. fj4Tl) and sums are reduced accordingly, then symmetric terms are 
collected into single sums by re-indexing the summations, obtaining 

^ j<k<l 



dVijki 

(43) 



The corresponding expression for the gradient of AGgb is similar but employs the Uij factors 
of Eq. (78) of reference '^^ rather than Wij. 

B.2 Component of derivative from Eq. (1391) 

Inserting Eq. (139|) in Eq. (63) of reference^ gives 

which is the same expression as that for the gradient of AG^av (see Appendix A of reference^) 
with the replacement 

Ik -r-WkPk 
An 
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The corresponding expression for the gradient of AGgb follows from the substitution: 



Ik -j-UkPk 
An 

B.3 Derivatives of HB Correction Energy 

From Eq. fl35|) we have 



Y.hsS'iws)^. (44) 



Fj g or 

Inserting Eqs. ([33]) and ([31D in Eq. (jS]) gives 



where 



dVsjk... _ ( dVsjk...\ dvsdVsjk 



where the first term on the r.h.s. represents the derivative of the overlap volume with respect 
to the position of atom i keeping the position of the water site s fixed, and the second term 
reflects the change of overlap volume due to a variation of the position of the water site 
caused by a shift in position of atom i. The latter term is non-zero only if i is one of the 
parent atoms of the water site. 
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